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The tremendous development of cold atom physics has opened up a whole new opportunity to 
study novel states of matter which are not easily accessible in solid state systems. Here we propose 
to realize the /-wave pairing superfluidity of spinless fermions in the pj;,y-orbital bands of the 
two dimensional honeycomb optical lattices. The non-trivial orbital band structure rather than 
strong correlation effects gives rise to the unconventional pairing with the nodal lines of the /-wave 
symmetry. With a confining harmonic trap, zero energy Andreev bound states appear around the 
circular boundary with a six-fold symmetry. The experimental realization and detection of this 
novel pairing state are feasible. 

PACS numbers: 



The study of unconventional Cooper pairing states 
has been a major subject in condensed matter physics 
for decades. In addition to the isotropic s-wave pair- 
ing, many unconventional pairing states have been 
identified-^. For example, different types of p- wave pair- 
ing states were found in the superfluid ^He systems^i^ in- 
cluding the anisotropic chiral ^-phase and the isotropic 
-B-phase. Evidence of the p-wave paring was also found 
in the ruthenate compound of Sr2Ru04^i^. The d-wave 
pairing states are most convincingly proved in the high- 
Tc cuprates with phase sensitive measurements including 
Josephson tunneling junction-'^ and zero energy Andreev 
bound states^i^. Many heavy fermion compounds exhibit 
evidence of unconventional Cooper pairing such as UPts, 
UBei3, and CeCoIus with nodal points or line s^^d"'^'"'^^ . 
However, phase sensitive measurements are still lacking. 
These unconventional pairing states are driven by strong 
correlation effects, which brings significant difficulties for 
theoretical analysis and prediction. It would be great to 
find novel unconventional pairing mechanisms easier to 
handle. 

On the other hand, cold atom physics has recently 
become an emerging frontier for condensed matter 
physic o-'^^d^'^^d^ . In particular, the s-wave pairing super- 
fluidity of fermions through the Feshbach resonances has 
become a major research focus. The Bose-Einstein con- 
densation (EEC) and Bardeen-Cooper-Schrieffer (BCS) 
crossover has been extensively investigate d^'^'d^dQi^Pi^i . 
Naturally, searching for unconventional Cooper pairing 
states in cold atom systems is expected to stimulate 
more exciting physics. For example, the p-wave pairing 
states have been proposed by using the p-wave Feshbach 
resonancea22i22ii^, which, however, suffer from a draw- 
back of heavy particle loss. The unconventional Cooper 
pairing states with cold atoms have not been realized 

In this article, we propose a novel /-wave pairing state 
of spinless fermions in the cold atom optical lattices. This 
unconventional pairing arises from the non-trivial band 
structure of the Pa,_y-orbital bands in the honeycomb opti- 
cal lattices combined with a conventional attractive inter- 



action. The internal orbital configurations of the Bloch 
wave band eigenstates vary with the crystal momenta, 
resulting in an /-wave angular form factor for the pair- 
ing order parameters. Along three high symmetry lines 
in the Brillouin zone whose directions differ by 120° from 
each other, the intra-band pairing order parameters are 
exactly suppressed to zero. The unconventional nature of 
this pairing exhibits in the appearance of the zero-energy 
Andreev bound states at the circular boundary with im- 
posing a confining trap. Since no strong correlation ef- 
fects are involved, our analysis below is well-controllable. 
This is a novel state of matter, which to our knowledge 
has not been unambiguously identified before neither in 
condensed matter nor in cold atom systems, thus this re- 
sult will greatly enrich the study of unconventional pair- 
ing states. 

The honeycomb optical lattice was experimentally con- 
structed quite some years ago by using three coplanar 
blue detuned laser beams whose wavevectors ki (i = 
1,2,3) differ by 120° from each other^^^. After the low- 
est s-orbital band is fulfilled, the next active ones are 
Px,j/-orbital bands lying in the hexagonal plane. Dif- 
ferent from the situation in graphene whose ^y-orhital 
bands strongly hybridize with the s-orbital bands and 
are pushed away from the Fermi surface, the hybridiza- 
tion between the Px,y and s-orbitals in the optical hon- 
eycomb lattice is negligible. The p^-orbital band can be 
pushed to higher energies and thus unoccupied by im- 
posing strong confinement along the z-axis. The P'x,y- 
orbital bands exhibit the characteristic features of or- 
bital physics, i.e., spatial anisotropy and orbital degen- 
eracy. This provides a new perspective of the honey- 
comb lattice, and is complementary to the research fo- 
cusing on the pz-hand system of graphene which is not 
orbitally active. The novel physics which does not ap- 
pear in graphene includes the flat band structure'^'', the 
consequential non-perturbative strong correlation effects 
(e.g. Wigner crystal^! and ferromagnetism^^) , frustra- 
tions in orbital exchange^, and the quantum anomalous 
Hall eSeciM. 

We employ the p^ .y-orhital band Hamiltonian studied 
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in Ref4a^>^^ as 

-^0=^11 X! {Pf-.i-Pr+ae^.j + h.C.} - ^ ^ Up, {^) 
feA, 1=1,2,3 reA<SB 

where ei.2,3 are unit vectors from one site in sublattice 
A to its three neighboring sites in sublattice B defined 
as ei,2 = i^ej; + ^Cy, 63 = -iy-, pi = {p^Cx + Pyiy) ■ ei 
is the p-orbital projected onto the bond along the direc- 
tion of i-i and only two of them are linearly independent; 
np = np^x + np^y is the particle number at site r; the 
(7-bonding describes the hopping between p-orbitals 
on neighboring sites parallel to the bond direction and is 
rescaled to 1 below; a is the nearest neighbor distance. t|| 
is positive due to the odd parity of p-orbitals. Eq. [1] ne- 
glects the TT-bonding ij^ describing the hopping between 
p-orbitals perpendicular to the bond direction. Due to 
the high spatial anisotropy of the p-orbitals, t_\_/t\\ <C 1. 
To confirm this, we have fitted tj_ and t\\ from a realis- 
tic band structure calculation for the sinusoidal optical 
potential V(r) = Vq X)i<i<j<3 cos{(fci - kj) ■ r\. With 
a moderate potential depth of Vq/Eh = 12, t^/t\\ is al- 
ready driven to 1% with t\\ « 0.?,7bER where Er = 
is the recoil energy and M is the atom mass. Further 
increasing Vq/Eji decreases tj^/ty even more. 

Eq. [T]has a chiral symmetry, i.e., under the transfor- 
mation of pp^x(y) ~* ~Pr.x(y) Only for sitcs in sublattice 
A but not for sites in sublattice B, Hq — > ~Hq. Thus its 
four bands have a symmetric spectra respect to the zero 
energy as 



^1,4 = T2*|h ^2,3 = /3 + 2^cosfc- (e,-e,) .(2) 

The Brillouin zone (BZ) is a regular hexagon with the 
edge length Two middle dispersive bands 2 

and 3 have two non-equivalent Dirac points K(K') — 
(±^^,0) with a band width of The bottom 

and top bands are flat. We define the four com- 
ponent annihilation operators in momentum space as 
(kik) = [pA,x{k),pA,yik),pB,x{k),pB,y{k)]. In this basis, 
the eigen-operator ipmik) for band m can be diagonal- 
ized as if^mik) = 4>n{k)Unm{k) where U{k) is a 4 x 4 
unitary matrix. The phase convention for band eigen- 
vectors, i.e., each column of U{k), is conveniently chosen 
as R:^il^m{k)R2^ = sgn{m)2pjn{k') with sgn(m) = — for 
TO = 1, 2 and sgn(TO) = -|- for m = 3, 4, where the sym- 
metry operation R:^ is the 60° rotation around a center 

of the hexagonal plaquette and k' = Ril k. The analytical 

form of Unm{k) is given in the Appendix A. 

The orbital configurations of band eigenvectors have 
interesting patterns as depicted in Fig. [2 which arise 
from the lattice Dg/i symmetry. Only the lowest disper- 
sive band (to = 2) is plotted as an example. The lowest 
fiat band (to = 1) has the same symmetry structure ex- 
cept having different orbital polarization directions. The 




FIG. 1; Orbital configurations of the lower dispersive band 
2 at high symmetry points and lines in the BZ. The time- 
reversal partners i/)2(±fc) only involve the same polar orbital 
if k is along the three middle lines, and take the complex 
orthogonal orbitals of ± ipy with opposite chiralities at K 
and K' . to = e'T . 



remaining two can be obtained by performing the opera- 
tion of the chiral symmetry. Six lines in the BZ are with 
the reflection symmetry, i.e., three passing the middle 
points of the opposite edges and three passing the op- 
posite vertices of K and K' . The eigenvectors of each 
band along these lines should be either even or odd re- 
spect to the corresponding reflection. For example, for 
the reflection respect to the j/-axis, the orbitals trans- 
form as Pa(B),x -PA(B).x,pA(B).,y ^ PA{B),y, thuS the 

eigen-operators il>m{k) with k \\ y must be either purely 
Py (even) or p^ (odd). For the other two middle lines, 
the corresponding eigenvectors are obtained by perform- 
ing the ±120° rotation. For the time-reversal partners 
i'm{k) and ipm{—k) along these fines, they only take the 
same real polar orbitals. On the other hand, the reflec- 
tion respect to the x-axis gives rise to the transforma- 
tion PA(B),x Pb(A),x and PA(B),y -PBiA),y Fur- 
thermore, K and K' have three-fold rotational symme- 
try. Combining two facts together, ij^miK) should be of 
Px + iPy for sites in one of the sublattices and Px — ipy 
for sites in the other sublattice. Its time-reversal partner 
ipmiK') has the opposite chiralities in both sublattices 
respect to those at K. In other words, the orbital con- 
figurations of ipmiki) are linear-polarized at the middle 
points of the BZ edge, changes to circularly-polarized at 
the vertices, and elliptically polarized in between. 

Next we introduce the on-site attractive interaction 
term between spinless fermions as 

Elint ^ U^^TlpxI^p^y, (3) 



where U is positive. We perform the mean-field decom- 



3 




FIG. 2: A) The /-wave pairing form factor F{k) for the intra- 
band pairing in momentum space. B) The /-wave pairing 
pattern in real space with Aa = —As = A. 
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J2 E Kmik)Mk)4'm{-k) + h.C, (4) 



FIG. 3: Gap and superfluid density for the /-wave state at 
l7/t|| = 3 and 0.5 < n < 1.5. A) The pairing amplitude 
A,4 = — A_g = A; B) and C) the lowest energy branch of 
Bogoliubov excitations ai n — 0.96 and n = 0.79. For a 
better visual effect, the negative eigenvalues are shown. D) 
The superfluid density ps v.s. n. 



where the pairing order parameters in the A and B- 
sublattices are self-consistently defined as A 



MB) 



\G)] {G\..\G) means the average over 

the pairing ground state; the summation of k only cov- 
ers half of the BZ; the multi-band pairing order param- 
eters in momentum space has a 4 x 4 matrix struc- 
ture as Anm{k) = U{G\'ipn{k)il'm{—k)\G). Under the 
rotation of it transforms as A„m(fc)_R^^ = 

sgn(n) sgn(m)A„„(/c'). 

The intra-band gap functions A„„ can be calculated 
analytically as A„„(A;) = i(— )"i(A^ — AB)F{k) with 
the /-wave form factor of 



F(k) = ^ sin ^^kx 

V3iVo(fc) 2 



cos — — kx 
2 



cos -ky 



, (5) 



where No{k) = |{3 - I]i<i<j<3 cos A: • {e, - ij)}. Aa 
can be fixed positive, and A^ — lA^le*'^^ with a relative 
phase A6. The optimal A6 takes the value of tt, i.e., 
Aa — —Ab, to maximize the intra-band pairings. We 
have confirmed this in the explicit self-consistent mean- 
field solution for Eq. [1] and Eq. H) Furthermore, the 
non-vanishing 7r-bonding tj_ term can further stabilize 
this solution as a result of the odd parity of 7r-orbitals. 

The configuration of Aa = —Ab = A exhibits the 
/-wave symmetry, i.e., the rotation of R:^ is equivalent 
to flipping the sign of A as depicted in Fig. [5] A. In 
particular, the diagonal terms satisfy i?| A„„(A:)i?^^ = 

— A„„(fc) as exhibit in F{k) with three nodal lines of 
k^ = 0, ky = ±fc^/V3. These are the same middle lines 
marked in Fig. [1] along which the time reversal partners 
Tpni^k) have the same polar orbital configurations. The 



pairing amplitude vanishes along these lines because in- 
teraction only exists between two orthogonal orbitals. On 
the other hand, maximal pairings occur between K and 
K' with opposite signs whose orbital configurations are 
orthogonal to each other. Unlike other examples of un- 
conventional pairing in condensed matter systems, this /- 
wave pairing structure mainly arises from the non-trivial 
orbital configuration of band structures but with conven- 
tional interactions. 

We next study the pairing strength A^ = — As = A in 
the weak coupling regime with U /t\\ = 3 in which the va- 
lidity of the self-consistent mean-field theory is justified. 
The lattice chiral symmetry ensures that these results 
are symmetric respect to the filling n = \. The case of 

< n < 0.5 corresponds to the filling in the bottom flat 
band in which each band eigenstate can be constructed 
as localized within a single hexagon plaquette presented 
(see Refi^i^). The attractive interaction drives all the 
plaquette states to touch each other leading to phase sep- 
aration until the flat band is fulfilled, which will be dis- 
cussed in a later publication. In Fig. [3] A, we plot A 
for 0.5 < n < 1.5 which corresponds to fill in the disper- 
sive bands. The maximal A appears dX n = 0.5 because 
all the particles in the flat band participate pairing. As 
n increases, A drops because pairing is less efficient in 
the dispersive bands. A reaches the minimum at n = 1 
where the Fermi energy touches the Dirac cones. 

Due to the multi-band structure, this /-wave pairing 
state remains fully gapped for general values of f//^|| and 
n. We plot the branch of the lowest energy Bogoliubov 
excitations for U/t\\ = 3 and the filling 0.5 < n < 1 
where the Fermi energy lies in band 2. For n close to 

1 in Fig. [3] B (n = 0.96), the Fermi surfaces form two 
disconnected pockets around the K and K' away from 
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FIG. 4: (a)The zero energy Andreev bound states appear on 
the boundary perpendicular to the anti-nodal directions, (b) 
The zero energy LDOS in real space with the maximal values 
a,t 9 — mr/S is a signature of the /-wave pairing symmetry. 



the nodal lines of A22(fc), thus the spectra are gapped. 
As lowering n in Fig. [3]C (n = 0.79), the Fermi surface 
becomes connected and intersects with the nodal lines 
of A22. At these intersections, the system in general re- 
mains gapped because of the non-zero inter-band pairing 

Ai2(fc). 

The superfluid density ps is plotted in Fig. [3] D for 
0.5 < n < 1.5 and C//i|| = 3, which is defined as 

Pt = = lim „ „ „ D 

^ 2m* 59^0 2 dSe^ ^ ' 

where 59 is the phase twist across the system boundaries. 
The behavior of ps is very different from that of the pair- 
ing gap, which is mostly determined by the states in the 
dispersive band. Close to n = 0.5, Ps is small in spite 
of the large value of A because of the immobility of the 
flat band states, ps reaches the maximal value about 
0.02t|| because the filling is close to the van Hove singu- 
larity of the density of states, and then drops as close to 
the Dirac points at n = 1. In two dimensional systems, 
the superfluidity develops below the Kosterlitz-Thouless 
(K-T) transition temperature Tkt ~ ^Ps- 

One of the most convincing proofs of the uncon- 
ventional pairing is the existence of the zero energy 
Andreev bound states at boundaries because of their 
phase sensitivity^. Fig. |4] A depicts the situation of 
the boundary perpendicular to the antinodal lines of 
the intra-band pairing. The scattering of the Bogoli- 
ubov quasi-particles changes momentum from fci„ to kref 
along which the pairing parameters switch the sign as 
^nn(kin) — — A„„(fcre/), which gives rise to zero energy 
Andreev bound states. On the other hand, if the bound- 
ary is perpendicular to the nodal lines, no phase changes 
occur and thus Andreev bound states vanish. 

Naturally in the experimental systems with an overall 
confining trap, the circular edge boundary samples all 
the orientations. We performed a real space Bogoliubov- 



de-Gennes calculation as 

( H A \ / u^{i) Unit) \ , . 

where H — Hq + Vex', Vex = ^^r^ is the harmonic 
confining potential with parameters Hfl = 0.02i|| and 

Iq = y/ ^ ~ 4.5a. The local density of states (LDOS) at 

energy E is defined as LD{f, E) = ^ J2n \un{r)\'^6{E — 
En) + \vnir)\'^d{E + En), where N is the total number 
of lattice sites. The numerical result of the zero energy 
LD{f,0) is depicted in real space in Fig. [3]B with the 
parameter values of = 3 and p = —0.18 which cor- 
responds to the filling in the center of the trap Ue = 0.97. 
Fig. [4]B shows that the LD{r, 0) is non-zero only at sites 
near the boundary, signaling the existence of the zero en- 
ergy Andreev bound states with six-fold symmetry. The 
LDOS maxima occur at angles 9 = n-| {n — 1 ^ 6) which 
are perpendicular to the antinodal directions shown in 
Fig. [HA, and the minima are located at 6* = (n + i)^. 

The experimental realization of this novel /-wave state 
is feasible. To enhance the attractive interaction be- 
tween spinless fermions, we propose to use atoms with 
large magnetic moments, such as ^^^Er with m = 7p,B 
on which laser cooling has been performed^. Compared 
to another possibility to use the p-wave Feshbach reso- 
nance, this method has the advantage to maintain the 
system stability. The interaction between two fermions 
in Px.y-ovhitals reads 

-U = y"dVidV2y(fi -f2){[l^p^(fi)Vp„(f2)]' 

- '^pArl)^Py{r2)1ppAr2)^Ppyirl)Y (8) 

where V'p^.p ^-re Wannier wavefunctions. The overall 
dipole interaction can be made attractive by polarizing 
the magnetic moments parallel to the hexagonal plane, 
which gives rise to V{fi—r2) = (1 — 3 cos^ 9) jr^ where 
r — |ri — ^2! and 9 is the angle between to and r\ — r-i. 
We take the laser wavelength A ~ GOOnm, and then the 
recoil energy E^ = 157nK. As shown in Appendix B, 
by choosing Vq/Eh = 30, the estimation shows that U w 
96nK and w 0.2^;^ = 31nK, and thus [//iy w 3 as 
chosen above. As shown in Fig. OD, the maximal Tkt ~ 
■|(0.02i||) « 1 ~ 2nK, which is within the experimentally 
accessible regime. 

A successful detection of the zero energy Andreev 
bound states will be a convincing proof to the /-wave 
pairing state. The radio-frequency (r/) spectroscopy has 
been an established tool to determine the pairing gap 
in cold atom system a^^i^^ . It also has a good spatial 
resolution"^^, which makes the direct imaging of the spa- 
tial distribution of the zero energy Andreev bound states 
feasible. The unique symmetry pattern of the zero energy 
Andreev bound states localized at the trap boundary can 
be revealed by identifying the locations of the zero energy 
spectrum determined from the spatially resolved rf spec- 
troscopy. 
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In summary, we have proposed the reahzation of a 
novel /-wave pairing state with spinless fermions which 
has not been identified in sohd state and cold atom sys- 
tems before. The key reason is the non-trivial p-orbital 
band structure of the honeycomb lattice rather than the 
strong correlation physics, which renders the above anal- 
ysis controllable. The Tkt is estimated to reach the 
order of InK within experimental accessibility. The rf 
spectroscopy detection of the six-fold symmetry pattern 
of the zero energy Andreev bound states along the cir- 
cular boundary will provide a phase sensitive test of the 
/-wave symmetry. 
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APPENDIX A: EIGENVECTORS OF THE BAND 
HAMILTONIAN AND THE PAIRING MATRIX 

In this section, we present the spectra of the band 
Hamiltonian Eq. [1] and the pairing matrix Amn ■ With 
the four-component spinor operator defined as 

= {PAx{k),PAy{k),PBx{k),PBv(k))'^, (Al) 

Eq. [1] becomes 

^0 = ti^J24>lik)Ho,mnik)Mk), (A2) 

k 

where the matrix kernel i?o,r?j.n(fc) takes the structure as 



where 



^^c'*'^! 



h.c. 



V3 ^gifc-ei 



v 








For each momentum fc, H()(k) is diagonalized as 

Ho,hn{k)Umn = EiUln- (A3) 

The band eigen-operators are expressed as 

V'm(fc) = 4>nik)Unrn{k). (A4) 

The unitary matrix U{k) reads 



Uik) 



( a*{k) e-'^b{k) e-'^h{k) a*{k) \ 

-b*{k) e-'^aik) e-''^a{k) -h*{k) 

a{k) e'"-^b*{k) -e'"-^b*{k) -a{k) 

V -b{k) e''-^a*{k) -e''-^a*{k) h(k) J 



(A5) 



a{k) 



/23(fc) -/3l(fc) . ^^^^^ /l2(fc) 



liN^{k) \/Noik) 



Mk) 
Noik) 



3 — E^ cos k ■ (Si — ij) 

l<i<j<3 



(A6) 



Using the eigenvector matrix U{k), the pairing matrix 
Amn can be spelled out straightforwardly as: 



{k) 



Ab 



For the solution of 
matrix as 



U2m{~k)Uin{k) - Ulm{~k)U2n{k) 
U4,ni-k)U3nik) - U 3m{~k)U in{k) . 

(A7) 

-Ab ~ A, we have the pairing 



A(fc) = 



/ Ai(fc) A2{k) Asik) \ 

A2(fc) -Ai(fc) -A3(fc) 

-A3(fc) -Ai(^) A2(fc) 

V Asik) A2(fc) Ai(fc) / 



,(A8) 



where the matrix elements read 
Ai(fc) 



16A . \/3, 
I — — — sm kx 

VSNoik) 2 



cos -r—kx — cos -ky 
2 2 ^ 



A2(fc) 
A3(fc) 



16A 

*3A^o(^) 

16A 
3A^o(fc) 



Ki{k) sin + K2{k) cos 



Ki (fc) cos Y - ^2 (fc) sin y 



(A9) 



and 



Xi(fc) = 4cos4^ + (4cos2^fc, 



7)cos2^ 
' 2 



V3, 



V3, 



K2{k) 



cos — COS -^fca; + 2 sin^ -^fcr 
2 2 2 

ky . \/3 J ky , 

sm —(cos — kr, — cos — 
2 ^ 2 2 



X [(2cos% + cos^fc^)^+sin2^fc^](A10) 



V3, 
2 



Both the band Hamiltonian Eq. [T] and the interaction 
Eq. [3] are invariant under the rotation R]l . Applying 
RjL to the band Hamiltonian matrix in momentum space 

Hoik), it transforms as 



-1 

3 



(All) 
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FIG. 5: F[y{h/l)] as a function of h/l. 



where Ril \s defined as 

3 



/ riL 



and 



Correspondingly, it can be shown directly that 
R^ipm{k) = sgn(m)'0m(fc') 



(A12) 



(A13) 



(A14) 



with sgn(TO) = —1 for TO = 1, 2 and +1 for to = 3, 4, and 
R^Amnik)R2^ = sgn(n)sgn(TO)A„„(fc'). 



APPENDIX B: CALCULATION OF THE 
ON-SITE HUBBARD INTERACTION WITH 
MAGNETIC DIPOLAR INTERACTION 

The optical potential around the center of each optical 
site can be approximated as 

V{r)^ju.\rl+rl) + '^c.lrl, (Bl) 

where we assume that the confinement in the z-axis is 
stronger so that ujz > OJ. The wavefunctions for the 
and py orbitals in this harmonic potential are 



'4'pAr) = C^Tj^expl - 
i^Pyir) = C^ry exp{ 



^2 1 2 2 



2P 



2li 



^2 I 2 2 
' X ^ y ' 



2P #}- 



where C = 21^H^^tt^^/'^ ,1 — ^hjmu, 1^ = ^JTiJrtuol, 
and Iz < I. 

With these wavefunctions, the on-site interaction can 



be evaluated with the direct and exchange terms: 

-U = j dPr,<fr2V{ri-r2)[[<PpAri)i^Py{r2)f 
-V'p. iri)ippy ir2)ipp^ (nj)V'p„ (f^i)} 

^ J d^rd^R V{r) (^!^kv±I^^ 

X Fi{r)F2{R), (B3) 

where we have introduced the relative coordinate r = 
ri — r2 and the center of mass coordinate R = (fi-|-r2)/2, 
and defined 



Fi{r) = exp{ 



2/2 



y _ ' z \ 



We propose to use fermionic atoms with large magnetic 
dipole moments to. By polarizing the magnetic moments 
with an external magnetic field, the anisotropic interac- 
tion reads 



V{r) 



fio — 3(to • f)]2 



(B5) 



where to = rn/m and f = f/r, respectively. We assume 
the polarization angle between to and the xy-plane is 6*, 
and perform the Gaussian integrals in Eq. IB3I The result 
is expressed analytically as 

^u = SiJd-^os^Wy), (B6) 

where 



(l+Z/)3/2[(3y+l)tan- 



(B7) 



and y is a function of l^/l as y{lz/l) — — iD- 

F[y{lz/l)] is a monotonic decreasing function of Iz as 
shown in Fig. [51 which reflects that the wavefunctions 
with a larger Iz have smaller overlaps with each other. 

The largest attractive interaction can be achieved by 
polarizing the magnetic moments in the xy-plane. If we 
use the fermionic atom of ^^^Er with m — 7^b and the 
laser beam with the wavelength A ~ 600nm to construct 
the honeycomb optical lattice, the recoil energy Ej^ — 
h'^/2m\'^ = 157nK and I « 93(-^)i/'*nm. Compared to 
the realistic band structure calculation with Vq /Er — 30 
in Ref.33, t\\ is fitted as ^ « Q.2Er = 31nK and I « 
40nm. With the choice of Iz/l — 0.2, U reaches reach 
96nK, thus U/t\\ sa 3 employed in this paper is achieved. 
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